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Abstract 

Four different kinds of laminar flows between two parallel plates are 
investigated using the Lattice Boltzmann Method (LBM) . The LBM accu- 
racy is estimated in two cases using numerical fits of the parabolic velocity 
profiles and the kinetic energy decay curves, respectively. The error rela- 
tive to the analytical kinematic viscosity values was found to be less than 
one percent in both cases. The LBM results for the unsteady development 
of the flow when one plate is brought suddenly at a constant velocity, are 
found in excellent agreement with the analytical solution. Because the 
classical Schlichting's approximate solution for the entrance-region flow 
is not valid for small Reynolds numbers, a Finite Element Method solu- 
tion was used in order to check the accuracy of the LBM results in this 
case. 



keywords: laminar flow, parallel plates, velocity profile, Lattice Boltz- 
mann. 



1 Introduction 

Since Frisch, Hasslacher and Pomcau jl], ||, have shown that particles mov- 
ing on a hexagonal lattice with very simple collisions rules on its nodes lead 
to the Navier-Stokes equation, the use of lattice gas models (LGM) to study 
hydrodynamics has received considerable interest. 
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The Lattice Boltzmann Method (LBM) ||, |, j| , a derivative of the lattice 
gas automaton method, uses real numbers instead of bits to represent particle 
distributions, being far less noisy than LGM. The LBM is a finite-difference 
technique for solving the kinetic equation in discrete space and discrete time; the 
Navier-Stokes equation is recovered in the long-wavelength and low-frequency 
limit. 

This paper presents a LBM study of four different laminar flows between two 
parallel plates. Both steady and unsteady flows are simulated; the numerical 
results are compared with analytical solutions or Finite Element Method (FEM) 
results. 



2 Lattice Boltzmann Method and Single Time 
Relaxation Approximation 

In the LBM, the evolution equation for the particle distribution functions rii(x, t) 
corresponding to the velocities c, , i = 0, 1, . . . 6 (cq = 0, | Cj \ = 1, j = 
1, 2, ... 6) defined on a two dimensional hexagonal lattice can be written as 
follows: 

rii(x + Cj, t + 1) = m(x, t) + Q.(rii(x, t)) (1) 

where fli = £l(rij(x,t)) is the local collision operator, depending only on the 
local particle distribution functions m. 

Chen et al. |Bj used the single time relaxation approximation 

n, = (2) 

where the equilibrium distribution function n^ q depends upon the local fluid 
variables and the lattice relaxation time r controls the evolution to the equilib- 
rium state. 

The total mass n and local momentum nu are conserved after collisions: 

i i 

nu = Ci7ii = ^ Ci n t 9 ( 4 ) 

i i 

The volumic density p for a hexagonal lattice is expressed in terms of n @: 

To derive the mass and momentum conservation equations, a Taylor expan- 
sion in time and space of eq.dJ) is performed in the long- wavelength and low- 
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frequency limit. Using the Chapman-Enskog procedure, the final results are: 



dn d 

m + d^ {nUa) = 



d(nu a ) d 
dt dx/3 



2t-1 d 
2 dxZ 



Ee 
Ci aCi flCi 







For the hexagonal lattice we have 

Ci = (Cii,C i2 ) 



n(i — 1) n(i — 1) 
cos , sm 



1 — do . I,-, . 2 2 
6 + g(Ci • u) + -(a ■ u) 

tiq 9 = n [do — u 2 ] 
where do S (0, 1) is a constant. This leads to: 



u 

7T 



1,2. ..6 
1,2. ..6 



E 



ra(l - d ) 



CipCi-fTli = -S a/3 + -rrc(V • m)(5 q/3 



(6) 
(7) 

(8) 

(9) 
(10) 

(11) 
(12) 



where S a p = § (§| 



is the time-rate-of-strain tensor. 

Dividing ([?]) with -\/3/2 and identifying the pressure as p = p(l — do)/2, the 
Navier-Stokes equation in the incompressible limit (V • u = 0) is obtained as: 



' V- [pl + puu pS =0 



from which we can identify the shear r) and kinematic v viscosities: 

2r - 1 , 2r - 1 

rj = p and f = 



S 



(13) 



(14) 



At the beginning of each run, the fluid was always considered to be at rest 
(u = 0). By choosing p = 2.1/(^/3/2) and d = 1/7, we have an uniform 
particle distribution m = 0.3, i = 0, 1, 2 ... 6 at each node. For each time step, 
the following successive operations are performed at every lattice node: 

• the local velocity is computed from (|j); 

• the equilibrium distribution functions n^ 9 (x, t) are computed from @ and 

Tol); 



• using (Q) and the propagation step is performed in order to obtain 
the rii(x + Ci, t + 1) distributions; 

• if a body force exist, the new m distributions are modified according to 
the induced particle momentum change during one time step. 
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Figure 1: The flow domain and the hexagonal lattice 




3 Flow between parallel rest plates 

Fig |l|. shows the flow domain between two parallel plates and the hexagonal 
two-dimensional lattice. The lattice is periodic in the flow direction, therefore 
the right nodes, connected with broken lines, coincide with the left ones. If Ni is 
the number of nodes along the x direction, the domain length is L — Ni because 
the unit length is the triangle side. The distance between two node rows in the 
y direction is V3/2. Because of the "bounce-back" rules imposed in order to 
satisfy the no-slip conditions on both planes 0, the zero- velocity boundaries 
are located at a distance -\/3/4 from the lowest and highest rows. If we denote 
with N w the number of node rows in the y direction, the distance between the 
two real plates is H = N w 

3.1 Steady unidirectional flow 

If the local velocity vector has the same direction everywhere and is independent 
of distance in the flow direction, the convective rate of change of velocity will 
vanish and the acceleration of a fluid element will be du/dt, with only one non- 
zero component. In addition, if we consider a steady flow, even this acceleration 
component will be zero. 
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The motion between two rest parallel plates can survive the viscous dissi- 
pation of energy only if there is a continuous energy supply to the fluid by a 
pressure gradient Vp, which must be also independent of x. When Vp is neg- 
ative, the pressure gradient represents an uniform body force per unit volume 
pf in the direction of the positive x-axis || . 

The x component of the equation of motion reduces to: 

3 11 

f + "W = (15) 

The solution which satisfies the boundary conditions u(0) = and u(H) = is: 

v 2 \H H 2 
with the maximum value 

f H 2 

u max = 0.125- (17) 

v 

Both the equations (fl5| ) and ( |l7j ) can be used in order to make a comparison 
between analytical (eq.|l4|) and numerical results. As an example, the kinematic 
viscosity v can be calculated by using the second derivative of u(y) from ( |l5| ) or 
the maximum velocity u max from ( |l7| ) . The first possibility was prefered in this 
paper because it takes into account the whole velocity profile shape. If we use 
the least square method to fit a parabola between (j/j, itj) points, i = 1 . . . N w , 
we will have for the second derivative 



dy 2 }2i{yi H ~ V 



I = - 2 ^7^— ^ (18) 



The body force per unit volume, pf, must equalize the decay of the momentum 
pu during one time step 

Taking into account that St = 1 for LBM, and 

where 5n is the value to be added to n±, respectively subtracted from at each 
step, according with (||) we obtain: 

/ = ^ (21) 

n 



(22) 



Finally, from (|15|),(|18[) and ( |2l|) the kinematic viscosity is: 

i = Sn YJsH^yf? 
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Figure 2: Numerical values for the kinematic viscosity v, obtained by fitting the 
parabolic profiles for 0.8 < r < 3; the solid line corresponds to the analytical 
formula (fill). 




The numerical results were obtained using a lattice with Ni — 100, N w = 101 
(H = 87.4686). For r > 0.8 we considered Sn — 10~ 6 , and the numerical values 
for v are represented in Fig. 0, together with the analytical line (|l4|). The 
relative error does not exceed 0.8%. 

A special attention was given for 0.5 < r < 0.8, the results being presented 
in Table [j]. It should be noticed that 8n must be decreased when r, respectively 
v 1 diminish in order to have only positive rii values, for each node and direction. 

Table |] also shows the number of steps needed to reach the steady state and 
the corresponding maximum velocity u max (see eq.|l7j). The relative difference 
between theoretical and analytical values, e„ = \v an — v num \ /v an , does not 
exceed one percent. 

Three velocity profiles are represented in Fig.||. A very good agreement exist 
between the theoretical parabola (solid curve) and the numerical results. 
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Table 1: Kinematic viscosities (y an from (|l4| ) and v num obtained by fitting 
parabolic velocity profiles) for 0.5 < r < 0.8. 
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3.2 Entrance— region steady flow 

If we consider wind-tunnel conditions at x — 0, i.e. an uniform velocity U in 
the inlet section, the velocity profile must become parabolic far downstream 

^) = 6C/ (|-^) ( 23 ) 

An analytical approximate description of such experimental situation has 
been performed by Schlichting jjj. At the beginning, i.e. at small distances 
from the inlet section, the boundary layers will grow in the same way as along 
a flat plate at zero incidence. The resulting velocity profile will consist of two 
boundary-layer profiles on the two walls joined in the center by a line of constant 
velocity. Since the flow rate must be the same for every section, the decrease 
in the rate of flow near the walls which is due to friction must be compensated 
by a corresponding increase near the axis. Thus the boundary layer is formed 
under the influence of an accelerated external flow, as distinct from the case of 
the flat plate. At larger distances from the inlet section the two boundary layers 
gradually merge into each other , and finally the velocity profile is transformed 
asymptotically into the parabolic distribution (|23|). Schlichting estimated the 
inlet length as Ie — 0.04 _ff Re, where Re denotes the Reynolds number referred 
to the width of the channel, U H/v. A comparison of results obtained with 
Lattice Gas Method to those computed by Schlichling was presented in pJDfl . 

The present LBM simulation was performed using a lattice with N w = 210, 
therefore H = 201 %/3/2 = 174.04. The average velocity between plates was 
chosen U = 0.01 and the kinematic viscosity v = 0.125 (for r = 1.0). It results 
a small Reynolds number Re = 13.925 and Ie = 97. In this case, the Schlichting 
solution is not valid because the boundary layer approximation cannot be used 
and it is expected that the Ie value becomes larger. To overcome this situation, 
we choose L = 200. 

In addition, in the boundary layer theory it is assumed that for incompress- 
ible flow the velocity profiles cannot have overshoots within the boundary layer 
but the results obtained with the Finite Element Method for entrance-region 
flow Jll| do not confirm this assumption. Fig. [| presents the velocity profile 
at four different distances from inlet section, obtained with LBM and FEM, 
respectively. 

In order to verify if the domain length was large enough, the dimensionless 
velocity, u/U, variation versus x/H is presented in Fig.|[ It can be seen that the 
curves reach asymptotically the values corresponding to the parabolic velocity 
profile (||). 
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Figure 4: Velocity profiles at different distances x/H = 0.0718 (a), 0.1436 (b), 
0.2873 (c), 0.5746 (d) from the inlet section (solid curves - LBM results, dots - 
FEM solution). 
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Figure 5: Velocity variation along the channel at different distances from the 
lower plate (solid lines - LBM results, dots - FEM solution). 
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3.3 Unsteady unidirectional flow 

Let us consider an unidirectional flow for which the body force vanishes at the 
initial moment t = 0. The equation of motion is: 

du d 2 u 

m= v W (24) 

with the boundary and initial conditions: 

u(0,t) = U(H,t) = 0, 
u{H/2,Q) = U 

The solution is: 

u(y, t) = U exp (-^ 2 ^j sin (^-|) (25) 
The free decay of the total fluid kinematic energy can be expressed as: 

(26) 



E c (t) 1 V 

and it can be used to calculate the kinematic viscosity from the slope of the 
In [E c (0) / E c (t)) <-> t line. Such lines are presented in FigJ|, the numerical 
results being obtained for a lattice with N w = 101 and iVj = 100. 

Table |^ presents the numerical and analytical values for v, the relative dif- 
ference being below 0.2%. 

4 Time development of steady flow between 
parallel plates in relative motion 

Suppose now that the fluid and the two parallel plates are initially at rest state. 
After the lower plate is suddenly brought to the steady velocity U in its own 
plane, while the upper one is still maintained at rest, the governing differential 
equation is (|24|) again, with the boundary and initial conditions 

u(Q,t) = U, u(H,t) = 0, for <>0 

u(y,0) = 0, for 0<y<H 

The velocity distribution is given by § : 

»<«•" ~" I 1 - 1) tE 7 ° p ir' V w) si " ("I) < 27 » 

i— 1 



The LBM simulation was made for a lattice with N w = 101 and Ni = 100. 
The distance between plates was now H = 101 V3/2 - V3/4. The term y/3/4 
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Figure 6: Kinetic energy decay. 




Table 2: Kinematic viscosities [y an from (jl4|) and v num obtained from the kinetic 
energy decay) for 0.6 < r < 3.0. 
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comes from the fact that we imposed the velocity U = 0.01 on the lower plate 
(there was no "bounce-back" condition there). 

The time evolution of the velocity profile is represented in Fig. fj], at t = 
300, 1200, 4800 and 10000, the kinematic viscosity being v = 0.125. 

An excellent concordance can be observed between the numerical and ana- 
lytical results. The linear velocity profile corresponding to the steady Couette 
flow with no pressure gradient it is naturally obtained for t — > oo. 



5 Conclusions 

A numerical approach to four kinds of flow between two parallel plates using 
the Lattice Boltzmann Method is presented. The number of lattice nodes was 
chosen in order achieve an acceptable CPU time on Alfa DEC computers. The 
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numerical and analytical values of the kinematic viscosity were used to estimate 
the accuracy of the LBM simulations for both steady and unsteady unidirec- 
tional flows, the relative differences being found below one percent. 

The zero velocity conditions on the rest channel walls were imposed by using 
"bounce-back" rules. The moving plate condition was achieved by forcing the 
the equilibrium distribution functions for the corresponding nodes in order to 
satisfy the imposed velocity. 

Since the Schlichting approximate analytical solution for the entrance-region 
flow is inadequate at low Reynolds numbers, the comparison of the LBM results 
with a Finite Element Method solution was found to be appropriate. 

Finally, we found an excellent agreement between the LBM computed time 
evolution of the velocity profiles and the analytical solutions of the Couette flow 
without pressure gradient. 
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